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ABSTRACT 


This thesis estimates the frequency response of a network where the only data 1s the 
output obtained from an Autoregressive-moving average (ARMA) model driven by a 
random input. 

Models of random processes and existing methods for solving ARMA models are 
examined. The estimation is performed iteratively by using the Yule-Walker Equations 
in three different methods for the AR part and the Cholesky factorization for the MA 
part. The AR parameters are estimated initially, then MA parameters are estimated 
assuming that the AR parameters have been compensated for. After the estimation of 
each parameter set, the original time series is filtered via the inverse of the last estimate 
of the transfer function of an AR model or MA model, allowing better and better esti- 
mation of each model's coefficients. The iteration refers to the procedure of removing 
the MA or AR part from the random process in an alternating fashion allowing the 
creation of an almost pure AR or MA process, respectively. As the iteration continues 
the estimates are improving. When the iteration reaches a point where the coefficients 


converge the last MA and AR model coefficients are retained as final estimates. 
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I. INTRODUCTION 


In this thesis, an iterative approach 1s presented to estimate the frequency response 
of a network where the only data 1s the output obtained from an Autoregressive-moving 
average (ARMA) model driven by a random input. The AR parameters are estimated 
initially, then the MA parameters are estimated assuming that the AR parameters have 
been compensated for.. 

To find the frequency response of an ARMA network two problems have to be ad- 
dressed. One is due to the shortness of the observed data (time limitation) and the re- 
sulting distortion of the estimated correlation function. The second problem is due to 
the nonlinear combination of the coefficients of the MA part as they appear in the esti- 
mate of the correlation function. 

To achieve this, we estimate coefficients of MA part and AR part iteratively. We 
use the Yule-Walker Method for AR coefficients estimation and the Cholesky decom- 
position for MA part estimation. We assume that the ARMA network output, 1.e., the 
observed signal is produced by white Gaussian noise driving the network. Three differ- 
ent methods are used to estimate the AR coefficients. To minimuze the effect of the 
MA part in the correlation function, correlation lags greater than the correlation length 
of the MA part are used in estimating the initial AR coefficients. The iteration refers 
to the procedure, which removes estimated MA or estimated AR contribution from the 
output of the ARMA model in an alternating fashion. This allows better and better es- 
timates of the AR and MA coefficients as the iteration continues. For the iterative AR 
coefficient estimation three methods can be used, denoted by method 1, method 2 and 
method 3. Method 1, 2 and 3 use correlation lags starting at p+ 1, 0, and 0 respectively, 
where the first two methods use a Square matrix inverse and the third method uses a 
Pseudo matrix inverse. 

Introduction to models of random processes is presented in Chapter II. Existing 
techniques to solve for the parameters of ARMA models are presented in Chapter II. 
Chapter IV includes simulation results. Simulation studies employed the Matlab pack- 


age on the IBM PC/AT. Conclusions and recommendations are presented in Chapter 
Mi. 


| 
I. MODELS OF RANDOM PROCESSES 


A. AR PROCESS MODEL } 

We say that x(n) 1s an autoregressive process of order p, or simply an AR(p) process | 

if it satisfies the difference equation, | 
X (1) ayn)? ee Ee) err ea (2.1) 


where a, a. ..., a, are constant coeflicients, and e(77) is a pure random process [Ref. 1]. 


Equation (2.1) may be written in the following form, 


f 


x(71) = — > x(t — k) + e(n) (272) | 


k=l | 


A realization of Liq.(2.2) is mlustrated in Pigure 1. 


AR process x(n) 





Figure 1. Autoregressive Model of Order p 


The system transfer function [1(z) between the input e(7) and the output x(n) for 


the AR model shown tn Figure 1 is 


Pa (2.3) 


H@)= = 
AZ) lt az tag t+. tayz? 





It is required that A(z) has all its roots within the unit circle in the z-plane which 
guarantees that H(z) is a stable and causal filter. 

The form of equation (2.3) illustrates that AR models have finite poles but no Ze- 
roes. Hence, this model is sometimes called an All-Pole model. Because an AR model 
has no poles outside or on the unit circle, it has the strict minimum-delay property and 
hence is always invertible. In general, minimum delay means that the transfer function 
must have no poles outside the unit circle, but can have poles on the unit circle. Strictly 
minimum delay means that the transfer function has no poles outside or on the unit 
eicle. 

The AR model ts also called an Infinite Impulse Response (IIR) filter. According 
to definition (2.2), output x(n) depends on past values of the output and on the present 
input. Because of this, it 1s also referred to as a Pure Feedback system. 


The power spectral density for AR models is given bv 


2 
O 


Parlf) = faa rr . (2.4) 
where 
p 
Aff) =]+ Nae 
k=] 


and T is the sampling interval [Ref 2]. 


B. MA PROCESS MODEL 
The sequence x(n) is said to be a moving average process of order q (denoted by 


MA(q)) if it satisfies the difference equation, 
x(71) = boe() + biel — 1) + .... + bg8(n — @) (>) 


where 6,, 5, ...., ,are coefficients, and e(”) is a pure random process [Ref. 1]. 


Equivalently, we may write, 


q 
x(n) =) yen — k) (2.6) 
k=0 
We may say that the output of an MA model depends only on present and past 
values of the input, Le., there is no feedback in an MA model. 


A realization of Eq. (2.6) is illustrated in Figure 2. 


MA process x(n) 





Figure 2. Mioving Average Nlodel of Order q 


The system transfer function 11(z) between the input e(77) and the output x(n) for 
the MA process 1s 


H(z) = B(z) = by + 6,27! + bgz* +. + B24 Oa) 


This transfer function has q finite zeroes, but no poles. Hence, the MA model ts 
also called an All-Zero model. MA models are invertible if and only if B(z) has no zeroes 
outside the unit circle, nor on the unit circle. 

MA models are also called Finite Impulse Response (FR) filters. 


The power spectral density for an MA model 1s given by 
ye a: | BY)? (2.8) 


where 





q 
BU) = by + > bye TT 
k—! 


and T is the sampling interval. 


C. ARMA PROCESS MODEL 
We say that x(n) is an autoregressive-moving average process of order (p,q) or 


simply an ARMA(p,q) process if it is satisfies the difference equation, 
erate) + 0d dx pl) One(z)-b bje(e— 1)+....+bje(n—p) (2.9) 


where, again, a, ...., a, Dp, ...., b, are coefficients and e(7) is a pure random process [Ref. 
la 


Equation (2.9) may be written as, 


P q ere) 
x(n) =— > age(n— k) +) byeln — K) =) Iyeln— b) (2.10) 
k=1 k=0 k=0 
The assumption }, = 1 can be made without anv loss of generality because the input 
e(n) can alwavs be scaled to account for any filter gain [Ref. 2]. 
A realization of Eq. (2.10) 1s illustrated in Figure 3. 


The system transfer function for the ARMA process 1s given by 


B(z) I+ biz + b,2~° ++ b,z4 





H(z) = e1) 


A(z) Ltaqz tag t...taz? 


where A(z) is the z-transform of the AR part and B(z) is the z-transform of the MA part. 

Both polynomials A(z) and B(z) are assumed to have all of their zeros within the 
unit circle of the z-plane to guarantee that H(z) is a stable minimum-phase invertible 
fier. 


The power spectral density for ARMA models is given by [Ref. 2] 


2 


BU) 
Aff) 


2 
ParmaY) = 4; 





(2.12) 








ARMA process x(n 





Figure 3.  Autoregressive-Moving Average Nlodel of Order (p,q) 


D. RELATIONSHIPS OF RANDOM PROCESSES 

The Wold decomposition theorem (Ref. 3} relates the AR, MA and ARMA models. 
It shows that, if the Power Spectral Density is purely continuous, any AR or ARMA 
process can be represented by a unique MA model of infinite order. 

Another important theorem which is stated by Kolmogorov [Ref. 4] says that any 


ARMA or MA process can be represented by an AR process of infinite order. 





To illustrate these theorems, we model [Ref. 5] an ARMA(I,1) process by an 
AR(co) or by an A(co) process. From equation (2.11), the system transfer function for 
the ARMA(I,1) process is 


lee bzs 


H(z) 
I + az" 


If we use AR(oo) process to represent ARMA(I,1) process, where 


l 
2 


Na 
Pelee ar Co ee 


and 


— feoace 
Cz) =1 +) qz7* = 


1+ 6,27 
=) 


By using synthetic division we find that 


C(z) = 1 + (a, — b,)z7' + (6? — a,b,)z7? + (6} — a, b{)27? + .. 


Hence inverse z-transform of az-” is ad(k — m) [Ref. 6] and 


; \ if k=0 
6(k) = 


0 else, 


the inverse z-transform of C(z) is 
c, = 5(k) + (a, — b,)d(k — 1) + (db; — @,b,)5(k — 2) + (bp — a,67)6(k — 3) +... 


Or 


if k=0 
a tb 6) a amit 1 


If we use a finite order AR(p) we should choose p to satisfv c,.,+0 or, equivalently 
b}x0. 
Therefore, a high-order AR model will be required when the zero of the ARMA 


process gets closer to the unit circle. 


In a similar way if we use an A{A(co) process to represent an ARMA(1,1) process, 
let 


H(z) = dp + ee + dz +... 


where 


= _ 146,27! 
D(z) = » az _ —_——— 
a ee ay2 


By using synthetic division as we did above, the inverse z-transform of D(z) will be 


: if k=0 
* ere| (Gra ce 


If we use a finite order MA(q) model, we should choose q to satisfy d,_,~0 or, 
equivalently a;x0. 
Therefore, a high-order MA model will be required when the pole of the ARMA 


process gets closer to the unit circle. 


FE. RELATIONSHIP OF AR, MA, AND ARMA PARAMETERS TO THE 
AUTOCORRELATION SEQUENCE 

In this section, we will present the relationship of the model parameters to the 
autocorrelation sequence [Ref. 2]. 


If we multiply Eq. (2.10) by x*(m — m) and take expectation, the result will be 


P q 


E{x(n)x*(n — m)} =—- > aE (x(n — k)x*(n—m)}+ > Eten — k)x*(n — m)} (2.13) 


k=] k=0 


where the superscript * is used to denote the complex conjugation. 


Equation (2.13) may be written as, 


P q 
rel) = — Y ayrgg(on — k) + Y Burl — b) 2.14) 
k=) k=0 


The cross correlation between the input and the output can be written as, 


ro(D = Efe(n + Dx*(n)} = Flan + nf + > Itet(n _ 0 || 


k=1 


rexll) = Yee) + Sick (+h) (2.15) 
k=1 
where A, = | by definition (Eq. 2.10). 
If we assume that the driving sequence is a white noise process of zero mean and 


variance o? then [Ref. 2] 


0 for 1>0 
r()=| of for 1=0 (2.16) 
op fon 7 =.0 


When we substitute Eq. (2.16) in Eq. (2.14), we get final relationship between the 


ARMA parameters and the autocorrelation sequence. 


Vy," ( —1) for m<0 
Pp q 
Z xe 
ma) = = > Airex (m— ky +02) yh ean for O<m<q (2.17) 
k=] p k=m 
- > apra(n — k) for m>q 
| k=1 


The relationship between the autocorrelation sequence and a pure autoregressive 


model may be written by setting g = 0 1n Eq. (2.17) 


Ya (m — k) for m>0 


: 
~ Litires Loli — 0 


ieee) for m<0 


ais) 


Yxx(m) = 


The relationship between the autocorrelation sequence and a pure moving average 


model may be wnitten by setting p = 0 in Eq. (2.17) 


0 for aa 


q 
rxx(M) = Ho Eee; fOr. Oe, (2.19) 
k=m 
rxx"( —m) for m<0O 


In this case we should note that 


h,=b, for leksq 


10 


Il. TECHNIQUES TO SOLVE FOR THE PARAMETERS OF ARMA 
MODELS 


The estimation of the parameters of ARMA processes 1s a classical problem which 
is still being investigated by statisticians. Methods to find the parameters for purely AR 
processes are well known, but for ARMA processes some problems remain. 

There is a nonlinear relationship between the ARMA parameters and the 
autocorrelation of process x(n). The nonlinear equation (2.17) presents the difficulty of 
estimating the ARMA parameters, even when we know the autocorrelation sequence 
exactly. Techniques based on iterative maximum likelihood estimation (MLE) can be 
used to find the ARMA parameters. These techniques require complex computations 
and are not guaranteed to converge, or they may converge to the wrong solution. 
Therefore they are not practical for real time series. For AR parameters, techniques 
based on the least squares criterion lead to solutions of linear equations and hence re- 
duce the computational complexity. Unfortunately, the moving average parameters of 
an ARMA model cannot be found easily by solving a set of linear equations. The MA 
parameters are convolved with the impulse response coefficients h(k) which causes a 
nonlinear relationship between the autocorrelation sequence and the filter coefficients. 

In section III-A and III-B, we will discuss the methods of AR and MA parameter 
estimation, and in section III-C we will present an iterative approach to find the pa- 


HalMete4rs. 


A. AR PARAMETER ESTIMATION 
In this section, we present three widely used methods of extracting the model pa- 
rameters from a given block of measured data x(n). 
These methods are: 
1. The autocorrelation, or Yule-Walker method 
2. The covariance method 


3. Burg’s method 


All three methods of estimating AR parameters are based on least-squares minimization 
criteria obtained by replacing the ensemble averages by appropriate time averages [Ref. 
7]. 


1] 


The criteria for the optimal forward ( e;(7) ) and backward ( e;(m) ) predictors are 


obtained by minimizing 
EL et(n)>] and E[ e(n)° | 


where e>(m) and e;(m) are the result of filtering x, through the prediction-error filter as 


given by 
e,(n) = xX, + Ap X,_) + pty 5 + Any X 
p — n P) n=) P2 n—2 Coeeeeseeees 
Cy (1) = Xp_p + Ap Xp—pgy TA Xn—pgz Peveerereeeee iat 3, 


The autocorrelation method 1s the most obvious and straightforward one. Equation 
(2.18) may be evaluated for the p+1 lag indices O< m< p and put into the following 


matrix form 


rex(O) Nye —1) re —2) re —P) IP y 5 

rl) r (0) a) asl pe Oa 1) ay 0 

Fe(2) Px) rxx(0) rf —P+2) 1121 =] 0 (3.1) 
rlP) telP— 1) relp — 2) r(0) IL%1 | 0 

By noting that 
Fe —K) = ryx(A) 
Eq. (3.1) can be written as 

70) een) rxx(2) relP) I y »? 

rel) —rxx(0) r<(1) rep — 1) {fa a 

rex(2) Taal 1) rx(0) exlP — 2) “2/=/ 0 (3.2) 
| rofP) relP 1) rel? —2) rea) |L®L LO 


Tron quo} 


P 


2 . 
o.= > relia, where a=1 
i=0 


We can replace the autocorrelations +7,,(k) by the corresponding sample 


autocorrelations (biased estimate) computed from the given block of data, where 


N=—1—k 


ro(k) = + > ox Oren Os hk =p (3.3) 


n=0 


We use the biased autocorrelation estimator because the unbiased autocorrelation 
estimates may result in autocorrelation matrices that are not positive semidefinite, which 
means certain matrix equations have no solution. On the other hand, the autocorrelation 
matrices formed from the biased autocorrelation estimate will always be positive semi- 
definite [Ref. 2 ]. 

Bv solving the normal equations, we can obtain estimates of the model’s parameters 
{d,, d, a,...., 4, 67}. Equation (3.1) is known as the AR, Yule-Walker or Normal 
Equations . The autocorrelation matrix of this equation 1s both Toeplitz and Flermitian 


because 


*x 


sed —k)= ot) 


The solution of the Hermitian Toeplitz equations can be computed with the 
Levinson Algorithm. This algorithm 1s a lattice realization for linear prediction filters. 
We illustrate the prediction-error sequence in Figure 4. 


The p” prediction error 1s given by 
€ (7) = Xyt Ay Xp + Ap Xpug Fees. + a,x (3.4) 
Looking at the first two prediction errors 
€;(") =X, Fay1%,_; 


e,(1) =x, + {Xp} ate 199% n—2 


where (1, a,,), and (1, @,,, a.) represent the best predictors of orders p=1 and p= 2, re- 
spectively. The extra index is used to indicate the order of the predictor. 
In the autocorrelation method the ensemble average in minimization of the forward 


prediction error is replaced bv the least-squares time average criteria @ as follows 


Is 


a Dya) SaoseOreeseoooue 





Figure 4. —‘ Prediction-error Sequence 


N+p—1 
ge Sa 2 
SC |e, (n)| 
A= 
N+p—1 (3,5) 
2 
Se 
= aa Xp + L_"pk*n-k 
R= 
n=0 

Where it 1s assumed that the data x, 1, ...... , Xy_, are observed and x, = 0 for outsidesmine 


range O<a<N—1. The minimization of the tume-average criteria with respect to the 
real and imaginary part of the a, ’s will Icad exactly to the same set of Yule-Walker 
Equations (3.1). One way to solve these equations is via the Levinson recursion which 
iS an iterative technique that extracts the next order predictor from the previous one. 
The Levinson algorithm can be summarized as follows [Ref. 7]: 
|- Initialize the recursion at p=0, by setting 
Ag(z) = | and Eg = r53(0) = Elan] 


where A,(z) 1s the prediction-error filter and £, is the mean-squared predic- 
tion error at the zero stage. So, initially we have no prediction. At stage p, 
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the prediction-error filter will be 4,(z) and the mean-squared prediction er- 
ror will be &. 


2- Compute Y,., which is called the reflection or PARCOR coefficient: 
P 
> arp = l = f) 
(=) 
Y 41 = ; (3.6) 
> ayr(0 
i=) 
3- Recursively determine the (p + 1)" order prediction-error filter polynomial 
FApa;(z); 
Be = Aa eet Aa(z *) (3.7) 
4- Update the mean-squared prediction error: 
2 
etme Lb n4)) => (3.8) 
5- Continue the iteration until the final desired order is reached. 


If the process x(n) is AR(p), then iteration will continue up to order p. It will pro- 
pigethe AR coefficients a,,, a,,, .... ,a,, which are also the best prediction coefficients. 
If we continue iteration after p, all prediction coefficients of order higher than p will be 
elese to zero {Ref. 7). 

Although the autocorrelation method is the most obvious and efficient one, and the 
resulting prediction-error filter 1s guaranteed to be minimal phase, it suffers from the 
effect of prewindowing the data sequence x(n) by padding it with zeros to the left and 
to the right. This reduces the accuracy of the method, especially when we have short 
data records. 

In the Covariance method the ensemble average in minimization of the forward pre- 


diction error is replaced by the time average as follows 


N-] 
Pp 
- 1 - 
i ye ps a > Apr nak i) 
k=) 
n=p 


The only difference between this method and autocorrelation method 1s 1n the limits 
of the summation. In the covariance method we observe all the data points needed to 


compute @. Since we do not need the data outside the range O< n< N—1 we do not 
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need to set x, explicitely to zero. For that reason the covariance method seems to be 
more realistic. 

If we perform the minimization, we find the AR parameter estimates as the solution 
of equations (3.10) [Ref. 5]. 


Ca a) os ,2) Geet! iP) Qy el 0) 
ev (2.1) mera) C,,(2,p) |] % ee (2.0) 
- - re: a (3.10) 
ea, 1) Cups) Cyx(Psp) “p Gortn ©) 
where 
NAA 
ex(j.k) = ae 
SD 
n=p 


Note that c,,(j,k) is an estimate of r,,g—k). The c,,(j,k) uses the sum of only N-p 
lag products to estimate the autocorrelation function for each lag even more lags are 
available. In contrast in the estimation of r,,(0) the autocorrelation method uses all data 
point, while the covariance method uses only N-p data points in the summation. The 
minimization gives us a non-Toeplitz c,,{j,A) matrix. This implies that we can not use 
the Levinson algorithm to solve Eq. (3.10). The equations may be solved by using the 
Cholesky decomposition which will be computationally more expensive. The estimated 
poles using this method are not guaranteed to lie within the unit circle. As an example, 
consider a first order predictor (p= 1) and a length-three (N= 3) sequence [Ref. 7]. The 


time-average criteria will be 


Z 
i 1] a 2 l 2 2 
? =>) Ie (n)| => (1 + ay Xo)" + (x2 + ay %)) 
n=! 
When we minimize @ by differentiating it with respect to a,, and setting the deriv- 
ative to zero , we have 
(x1 + 4, Xo)xXo + (Xq + yx )x1 = 0 


X1Xq + XX, 
C= 

2 2 
Xo + xy 
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Note the denominator does not depend on the variable x, and if x, is large enough 
we might have magnitude of a, greater than one. In this case we do not get a 
minimum-phase prediction-error filter which is sometimes not desirable 1n practice. 


Bureg’s method estimates the reflection coefficients and then uses the Levinson 





recursion to estimate the AR parameters. Burg’s minimization criteria 1s to minimize 


both the forward prediction error e>() and the backward prediction error e>(m): 


N=1 
A | = 
OS N=. L Let a (n)° | 


The computational steps are summarized below [Ref. 7]: 


1- Initialize, bv setting ; 
es (n) = eo (n) = Xx, for O<n<N-1 (3.11) 
and 
N-1 
LILY 2 
Eo = N Xn 
n=0 
At stage p-1, the prediction-error filter will be A,_,(z) which 1s the Z trans- 
iomanor tMcecc(uciice (lew ecrennenes) dj) lhe mean-squared error 
Timber ea), andwee. 7) ean be calculated for p-l<ns .\—1. 
2- Compute the reflection coefficient: 
N=1 
fe _ 
2) ef ane _,(2 — 1) 
ne ie Babe 
a i 
ae (3.12) 
+ 7.2, = 2 
» efit cee )) 
n=p 
To guarantee the minimum-phase property, Y, must have a magnitude less 
than one. 
3- Compute 4,(z) using the Levinson recursion given by: 
I ] 0 
An 1 An—1,1 An—1 p—1 
a a a 
P,2 p—1,2 p—1p-2 
ap p-l Gp—1,p—I Gp—,1 
Ap,p 0 1 
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4- Compute e>(m) and e-(m) forp<nsN—1. 


ep (n) = e7_4(7) — ei a) (3.14) 
a aa (a — leeg@ 
3- Update the mean-squared error as follows: 
(lee (3.15) 
6- Continue the iteration until p equals the model order. 


The Burg’s method estimates the poles which are on or inside the unit circle. This 
is due to the property | Y,| <1. Therefore, care must be taken to deal with the situation 


when | Y,] = 1, as this causes the prediction filter to become non-minimum phase. 


B. MA PARAMETER ESTIMATION. 

The most obvious approach to estimate the MA parameters would be to solve the 
nonlinear equation (2.19) using the autocorrelation sequence. Solutions of Eq. (2.19) 
involve difficult spectral factorization techniques [Ref. 8]. 

There is another approach called Durbin’s method which 1s related Kolmogorov's 
theorem [Ref. 4] and is based on a high order AR approximation of the MA process. 


The AR process allows results using only linear operations. Let 


q 
Be) = 1+ > 27 
k=) 


represent the system transfer function of an MA(q) process, and 


A(z) =1+ > az 
k= 


represent the system transfer function of an AR(oo) process that is equivalent to the 


MA(q) process. Therefore, we have 


l 


32) = 7 @) 





or 


B(z)A,_(z) = 1 (3.16) 


The inverse Z-transform of the Eq. (3.16) is: 


: (3.17) 


q 
1 for m 
Am + > Pndinn = 0(m) = ‘0 for m RE 


n=] 


where a, = | and a, = 0 for k < 0. Therefore, the MA parameters can be determined from 
the infinite-order AR model by solving (3.17). 

In practice, one can calculate high-order AR(M) parameters, where > q. Based 
on these parameter estimates (1, @,(1), @,,(2), ....., @,(49) , an error in the MA part is 
Sommputed {Ref. 2]. 


eral) = Gag) +) Bydyglm — 1) (3.18) 


—" 


According to Eq. (3.17) the error should be zero for all m except for m=0. But in 
practice, the error will not be zero when using finite data, so MA parameter estimates 


are obtained by the minimization of the squared error power, given by 


2 
A e 1 


m=0 


This estimation procedure is an approximate maximum likelihood estimate (MLE). 
The approximate MLE procedure using Durbin’s method for MA parameters results in 


the following estimation [Ref 5]. 
peas ae. (3.20) 
where 


M-— [i-/| 


& | A A _ 
[ Rea ly = Say by An ans |i—j| for sow == | neers +g 


n=0 


M-l 


“A | N A 4 
Cheeks) bad for lie 2 ese cans MG, 


7=0 


Ih) 


Again the Levinson algorithm may be used to solve Eq. (3.20) for the b parameters. 
The estimated zeros of B(z) will be inside the unit circle by the minimum-phase property 
of the autocorrelation method. 

In summary, Durbin’s method first uses the data x, x,,..-, Xy_, to find a largewemeen 
AR(M) model using the autocorrelation method. Then using these AR parameter esti- 
mates (I, ay ayaa as tessa, (b,, b,, ake b,) is found. 

Another technique for estimating the MA parameters is to use the MA spectral 


estimator [Ref: 5] which is given by Eq. (3.21) based on sample correlation values. 


q 


Pus = ) Fes(k)e 2 (3.21) 


k=—q 


Since theoretically Eq. (3.21) 1s equal to g3| Bf) | 2, the MA parameters can be found 
by using the Spectral factorization theorem [Ref. 9]. This theorem shows that anv ra- 
tional power spectral density of a stationary signal x, can be factored into a minimum- 


phase form 
P,x(2) = 0; B(z) Bz") (3.22) 


C. ITERATIVE ARMA PARAMETER ESTIMATION 
Let us consider the modelling of ARMA transfer function as given in Figure 5. 
The computational steps are summarized below for the iteration assuming know- 


ledge of the model order: 


I- Form the sample correlation of the time series y(n) by using the biased 
autocorrelation estimator. 


2- Get the AR coefficients using the correlation lags > g (to minimize the MA 
influence). 

3- Inverse filter the original time series y(n) via the inverse of the AR filter (use 
the AR coefficients of step 2) to get x,(n) 

4. Form the sample correlation of the time series x,(7) by using the biased 
autocorrelation estimator. 

5- Get the MA coefficients using the sample correlation of x,(m) . 

6- Inverse filter the original time series v(n) via the inverse of the MA filter 


(use MA coefficients which are estimated in step 5) to get x,(”). We assume 
that the MA coefficients are of minimum phase (all roots are inside the unit 
Cine): 
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Modelling of ARNIA transfer function. 


Form the sample correlation of time series x,(7) by using the biased 
autocorrelation estimator. 


Get the AR coefficients using all lags or Jags > g depending on the method 
sed: 


Inverse filter the original v(n) series using the AR coeflicient estimates 
Which are obtained in step 8, to get x,(7). 


Get the sample correlation of x,(7) by using the biased autocorrelation es- 
timator. 


Compute the MA coefficients using the sample correlation of x,(72) . 


Compute the error in the AR and the MA parts as given below 
p q 


—1 —1 
errorj — l)= ) (a _ ra \ + iG _ b, li fore) — 2, 5,00 10(3.23) 


i=l i=) 
Where superscript J is used to denote the number of iteration. The upper 


count of j is experimentally chosen. Iteration continue up to 10 if the co- 
efficients do not converge. 


If error > A, then go to step 6, else exit the program using the last updates 
of the filter coefficients. If > 10 terminate with an error message. Note 7 
is a small experimentally chosen number. 


Z| 


The AR coefficients can be obtained via a pseudoinverse from the modified Yule- 
Walker equations [Ref. 2]. Equation (2.17) may be rewritten for the p lag indices 


g+Il<m<q+p,and put into matrix form 


id) Gia) rAd — Pawel nay TG ae 
r r — 2 a r 
at 1) — Pe, — ) ‘ = —— 2) (3.24) 
rodQ@t pth) re(gtpt2) reolg)  |L® req +P) 


We can compute the autocorrelation sequence for lags q-p+1 to q+p, therefore 
the AR parameters are found as the solution of Eq. (3.24) where the MA parts influence 
is minimal. The autocorrelation matrix 1s of Toeplitz form. 

The Moore-Penrose pseudoinverse A* of an m x n matrix A provides the minimum- 


norm least squares solution to Ax =b. The solution is given by 
Se 


where x 18 am x 1 vector that simultaneously minimizes the squared equation error, for 
a given x | vector b [Ref. 2]. If m=n and rank A 1s n (1.e., A 1s nonsingular), then the 
pseudoinverse becomes the square matrix inverse A*=A™. If m>n (1e., more 


equations than unknowns) and rank A 1s n, then 
A*=(A" A) A" and 


x=(A"A)'A"b 


The superscript H 1s used to denote the Hermitian transpose operation. This is the 
least squares solution for a set of overdetermined equations. 

In step 2 of the iteration, we will use correlation lags greater than q while in step 
8 either all lags or lags greater than q are used depending on the method. 

For MA coefficients estimation the Cholesky decomposition is used. If the matrix 
A 1s square and Hermitian, then the usual triangular factorization takes on the special 


form 


A=RR"® (3.25) 
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where R is a lower triangular matrix with nonzero real principal diagonal elements. This 
decomposition is called the Cholesky decomposition (Ref. 2]. 


For the MA part, the statistical autocorrelation is given by 
q—\k\ 


Rylk) = 02 Y Pngubn (3.26) 


n=0 


Let B be the lower triangular Toeplitz matrix where the matrix elements 6,, are given 
in terms of the impulse response of the filter B(z) [Ref. 7]: 


and let the autocorrelation matrix of x, be 
Sell (5 = Hells a, 


Then, the transpose matrix B? will have matrix elements 


and Eq. (3.26) can be written as 


pase) 
Ryglisf) = Ryxli-D= 2 YD OngtyPr 


n=0 


Ryx(i,/) = 0;(BB)y 
Thus, in matrix notation 
R,,= 0. BB" (3.27) 


which is related to Eq. (3.25) with the assumption that 


Therefore, an approximation of the MA parameters can be found from the Cholesky 


decomposition of the correlation matrix R,,. 


Pedy 


IV. SIMULATION RESULTS 


In this chapter, computer simulation results are presented to show how three dif- 
ferent estimation methods work on various ARMA models. Two ARMA(2,2) models 
which differ in pole-zero locations, an ARMA(2,3) model and an ARMA(3,4) model are 
used as test models. In addition two realizations for each ARMA process are utilized. 
Two hundred data points are used in computing the sample autocorrelation values. The 
three different AR estimation methods are explained below. 

Method 1 

The AR coefficients are obtained via a square matrix inverse using the modified 
Yule-Walker equations. Correlation lags greater than q are used to minimize the MA 
part influence at the first calculation. For the remainder of the iteration the same cor- 
relation lags are used to minimize potential influence from the MA part. A Cholesky 
factorization is used to find the MA part coefficients. 

Method 2 

The AR coefficients are obtained via a square matrix inverse using the modified 
Yule-Walker equations. Correlation lags greater than q are used to minimize the MA 
part influence at first calculation. For the remainder of the iteration the correlation lags 
Starting from zero are used assuming the MA part contribution has effectively been re- 
moved. The Cholesky factorization is used to find the MA part coefficients. 

Method 3 

The AR coefficients are obtained via a pseudoinverse instead of a square matrix 
inverse using the correlation function starting at the zero lag after first calculation. This 
allows the use of all important correlation lags. The Cholesky factorization is used to 
find the MA part coefficients. 

In addition, observation noise is added to one of the ARMA models and the re- 
sulting noisy sequence is processed via the three different methods. The observation 
noise is independent of the driving noise. It is white Gaussian noise with zero mean and 
unit variance. The signal to noise ratio (SNR) is about 15 dB. 

The computer program computes the differential errors in the AR and MA parts 
as given in Eq.(3.23), and then compares it with a small experimentally established value 


/ (i.e., 2 = 0.0001). If the error is less than / the program is terminated. If the error 1s 
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larger than / the program 1s reentered at step 6. Also, if the coefficients do not converge 
the program will stop after nine iterations. 

Comparisons of the true and estimated coefficient differences, of pole-zero lo- 
cations, of distances between the true and estimated pole-zero locations and of radial 
differences between the true and estimated pole-zero locations are presented to show 
how well the three methods work. Also the spectra of the ARMA models are plotted by 


using the true and estimated coefficients. 


A. THE ARMA(2,2) MODEL-A 


The pole-zero locations for this-model are illustrated in Figure 6. 





Figure 6. The ARMA(2,2) Model-A, Pole-zero Locations 


— 


. METHOD 1 
a. Noise Realization I 
The coefficients converge after two iterations. Tables I and 2 present the 


results. 


» 


Table 1. ARMA(2.2) MODEL-A, METHOD 1, 
COEFFICIENTS COMPARISON, 
REALIZATION 1. 


Ta | 180 | 18412 | 00a? 
[| 1.00 | 1.3053 | +0053 
[4,080 | 0.6839 [0.1561 













Table 2. ARMA(2,2) MODEL-A, METHOD 1, POLE-ZERO COMPAR- 
ISON, REALIZATION 1. 


Zeros 
0.9206 + 0.2250 










b. Notse Realization 2 
Using a different noise realization, the coefficients still converge after two 


iterations. Tables 3 and 4 present the results. 


Table 3. ARMA(2,2) MODEL-A, METHOD 1, 
COEFFICIENTS COMPARISON, 
REALIZATION 2 


a | 180 | 187-0072 
[| 1.00 | 110s3_| +0.1083 
4, | oso | osm | -02929 
[tf oso | 0.1502 | 0.698 
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Table 4. ARMA(2,2) MODEL-A, METHOD 1, POLE-ZERO COMPAR- 
ISON, REALIZATION 2 


Poles- True Estimated Distance Radial Diff. 
Zeros 










The Spectra using the true and the estimated network coefficients are 
plotted in Figure 7. 
2, METHOD 2 
a. Noise Realization I 
The coefficients converge after two iterations. Tables 5 and 6 present the 


results. 


Table 5. ARMA(2,2) MODEL-A, METHOD 2, 
COEFFICIENTS COMPARISON, 
REALIZATION 1. 


3, | oso | o.sas7 011533 








2 


Table 6. ARMA(2,2) MODEL-A, METHOD 2, POLE-ZERO COMPAR- 
ISON, REALIZATION 1. 


Poles- True Estimated Distance Radial Diff. 
Zeros 












b. Noise Realization 2 
Using a different noise realization, the coefficients still converge after two 


iterations. Tables 7 and 8 present the results. 


Table 7. ARMA(2,2) MODEL-A, METHOD 2, 
COEFFICIENTS COMPARISON, 
REALIZATION 2. 


| 1.00 | 1.0982 | +0.0983 








Table 8. ARMA(2,2) MODEL-A, METHOD 2, POLE-ZERO COMPAR- 
ISON, REALIZATION 2. 


Zeros 
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The spectra using the true and the estimated network coefficients are plot- 
fedein Figure 8. 
3. METHOD 3 
a. Noise Realization I 
The coefficients converge after two iterations. Tables 9 and 10 present the 


results. 


Table 9. ARMA(2,2) MODEL-A, METHOD 3, 
COEFFICIENTS COMPARISON, 
REALIZATION I. 


[| 1.00 | 1.3062 | 0.3062 
4,080 [0.6382 [0.1518 













Table 10. ARMA(2,2) MODEL-A, METHOD 3, POLE-ZERO COM- 
PARISON, REALIZATION 1. 


Zeros 














b. Noise Realization 2 
For a different noise realization, the coefficients still converge after one it- 


eration. Tables 11 and 12 present the results. 


2 


Table 11. ARMA(2.2) MODEL-A, METHOD 
3, COEFFICIENTS COMPAR- 
ISON, REALIZATION 2 


a | 180 | 8545 | 0.0545 
a 
a a 








Table 12. ARMA(2,2) MODEL-A, METHOD 3, POLE-ZERO COM- 
PARISON, REALIZATION 2 


Poles- True Estimated Distance Radial Diff. 
Zeros 


0.2365 + 0.5260 








The spectra using the true and the estimated network coefficients are plot- 


ted in Figure 9: 


three methods. 


terms of the coefficients and pole-zero locations. 


In summary, the estimated coefficients converge after two iterations in all 


clients converge after one iteration using method 3. 
In each method the AR part coefficients of the ARMA(2,2) model tend to 


be more accurate than the MA part coefficients. 


THE ARMA(2,2) MODEL-B 


The pole-zero locations for this model are illustrated in Figure 10. 
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For both noise realizations, the third method gives the best result in 


For the noise realization 2, the coeffi- 


ARMA(2,2) Model A, Method 1, Reafizatiou If 






—— Using trite coefficients 


*#¢* ‘After 2 iterations 





- tot ep ne Cmneems SOE etemtnn som: 






a 
oO 
oO 
AS, 
5 
ee 
G 
aD 
E 
> 
Sample Number 
ARMA(2,2) Model A, Method 1, Realization 2 
0 t 
— Using true coefficients 

*##* After 2 iterations | 

—10]- js Lae | 

—_ 

| | 

~20 | 

a : 

B 72° | 

42 | : 

: 

2 : 

f —40}- | 

6 | 

= 

| 


-50 


=60 


0 10 20 30 40 50 60 70 


Sample Number 


Figure 7. The Spectra of ARMA(2,2) Model-A, Method ! 
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ARMA(2,2) Model A, Method 3, Reafization 1! 
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Figure 9. The Spectra of ARMIA(2,2) Model-A, Method 3 
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Figure 10. The ARMA(2,2) Model-B, Pole-zero Locations 


1 METHOD 1 
a. Noise Realization 1 
The coefficients converge after three iterations. Tables 13 and 14 present 


the results. 


Table 13. ARNIA(2,2) MODEL-B, METHOD 
1, COEFFICIENTS CONMPAR- 
ISON, REALIZATION 1. 






0.3125 -0.1317 -0.4442 
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Table 14. ARMA(2,2) MODEL-B, METHOD 1, POLE-ZERO COMPAR- 
ISON, REALIZATION I. 


Poles- True Estimated Distance Radial Diff. 
Zeros 






b. Noise Realization 2 
For a different noise realization, the coefficients converge after two iter- 


ations. Tables 15 and 16 present the results. 


Table 15. ARMA(2.2) MODEL-B, METHOD 
1, COEFFICIENTS COMPAR- 
ISON, REALIZATION 2 













Table 16. ARMA(2,2) MODEL-B, METHOD 1, POLE-ZERO COMPAR- 
ISON, REALIZATION 2 


Poles- True Estimated Distance Radial Diff. 
Zeros 


0.255 











fad | rd 


oD 


The spectra using the true and the estimated network coefficients are plot- 
ted in Figures ale 
2, METHOD 2 
a. Noise Realization 1 
The coefficients converge after six iterations. Tables 17 and 18 present the 


results. 


Table 17. ARMA(2,2) MODEL-B, METHOD 
2, COEFFICIENTS COMPAR- 
ISON, REALIZATION 1. 


£00969 
0.5013 | -000187 


| 4.00 | 1orma | +0074 
4, P0950] _0.7s2 [0.0268 
03128 





Table 18. ARMA(2,2) MODEL-B, METHOD 2, POLE-ZERO COMPAR- 
ISON, REALIZATION 1. 


Poles- True Estimated Distance Radial Diff. 
Zeros 


0.6 + 0.4) | 0.5766 + 0.4110} 0.0312 


0.6 - 0.44 0.5766 - 0.4110} 0.0312 
02336 + 02535 | 02473 
0.25 = 0.5) ||ROaearer see 0.2793 





b. Notse Realization 2 
For a different noise realization, the coefficients converge after five iter- 


ations. Tables 19 and 20 present the results. 
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Table 19. ARMA(2.2) MODEL-B, METHOD 
2, COEFFICIENTS COMPAR- 
ISON, REALIZATION 2 


er 
| 1.00 | 10610 | +0.0610 
4, | 050 | 0.1301 | 03609 








Table 20. ARMA(2,2) MODEL-B, METHOD 2, POLE-ZERO COMPAR- 
ISON, REALIZATION 2 


Zeros 
0.7957 + 0.3477} 








The spectra using the true and the estimated network coefficients are plot- 
eden figure | 2. 
3. METHOD 3 
a. Noise Realization 1 
The coefficients converge after four iterations. Tables 2] and 22 present 


the results. 
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Table 21. ARMA(2.2) MODEL-B, METHOD 
3, COEFFICIENTS COMPAR- 
ISON, REALIZATION 1. 


Ta | 120) -Lisis | + 00182 
am | 1.00 [1.0030 | +0.0030 








Table 22. ARMA(2,2) MODEL-B, METHOD 3, POLE-ZERO COMPAR- 
ISON, REALIZATION I. 


Zeros 






b. Notse Realization 2 
For a different noise realization. the coefficients converge after two iter- 


ations. Tables 23 and 24 present the results. 


Table 23. ARMA(2,2) MODEL-B, METHOD 
3, COEFFICIENTS COMPAR- 
ISON, REALIZATION 2. 


[| 080 | 0use7 | 03173 








38 


Table 24. ARMA(2.2) MODEL-B, METHOD 3, POLE-ZERO COMPAR- 
ISON, REALIZATION 2. 


Zeros 






The spectra using the true and the estimated network coefficients are plot- 
feain Figure 13. 

In summary, under each noise realization the method 3 gives the best result 
in terms of the coefficients and pole-zero locations. 

For each method, the spectrum of the ARMA(2,2) model using estimated 
coefficients indicates reasonably accurate pole locations but the zero locations of the 
spectrum do not follow the original ones. This means that the AR part coefficients of 


the ARMA model tend to be more accurate than the MA part coefficients. 


C. THE ARMA(3,4) MODEL 
The pole-zero locations for this model are illustrated in Figure 14. 
1 METHOD 1 
a. Noise Realization 1 
The coefficients do not converge in nine iterations. Because of oscillations 
about two values, the average of the two values is used as an estimate of the coefficients 


in the comparison. Tables 25 and 26 present the results. 
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Figure 14. The ARNIA(3,4) Model Pole-zero Locations 


Table 25. ARNIA(3,4) MODEL, METHOD 1, 
COEFFICIENTS COMPARISON, 
REALIZATION 1. 


Estimated Difference 

+0.0725 

05158 
| 





18 le 

E0335 
0.2509 
+ 2.3059 
1.3473 0.9473 
1.6482 -2.1282 
2.0981 
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Table 26. ARMA(3,4) MODEL, METHOD 1, POLE-ZERO COMPAR- 
ISON, REALIZATION I. 


Poles- True Estimated Distance Radial Diff. 
Zeros 









b. Noise Realization 2 
Using a different noise realization, the coefficients still do not converge in 
nine iterations. Because of oscillations about two values, the average of the two values 
is used as an estimate of the coefficients in the comparison. Tables 27 and 28 present 


the results. 


Table 27. ARMA(3,4) MODEL, METHOD 1, 
COEFFICIENTS COMPARISON, 
REALIZATION 2. 


re 
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Table 28. ARMA(3.4) MODEL, METHOD 1, POLE-ZERO COMPAR- 
ISON, REALIZATION 2. 


Poles- True Estimated Distance Radial Diff. 
Zeros 












The spectra using the true and the estimated network coefficients are plot- 
foam Figure 15. 
2, METHOD 2 
a. Notse Realization 1 
The coefficients converge after nine iterations. Tables 29 and 30 present the 


results. 


Table 29. ARMA(3,4) MODEL, METHOD 2, 
COEFFICIENTS 
COMPARISON,REALIZATION 1. 


[2 | 040 | 02668 | 0.1336 
oat [0.0262 [0.5002 
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Table 30. ARMA(3,4) MODEL, METHOD 2, POLE-ZERO COMPAR- 
ISON, REALIZATION 1. 


Poles- True Estimated Distance Radial Diff. 
Zeros 












b. Noise Realization 2 
Using the second noise realization, the coefficients converge after nine it- 


erations. Tables 31 and 32 present the results. 


Table 31. ARMA(3.4) MODEL, METHOD 2, 
COEFFICIENTS COMPARISON, 
REALIZATION 2. 


-1,6668 -0.0668 


0.7235 
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Table 32. ARMA(3,4) MODEL, METHOD 2, POLE-ZERO COMPAR- 
ISON, REALIZATION 2. 


Zeros 


0.6271 - 0.1572 





The spectra using the true and the estimated network coefficients are plot- 
ted in Figure 16. 
3. METHOD 3 
a. Notse Realization 1 
The coefficients converge after six iterations. Tables 33 and 34 present the 


results. 


Table 33. ARMA(3,4) MODEL, METHOD 3, 
COEFFICIENTS COMPARISON, 
REALIZATION 1. 


[00 fnims3_| +0733 
[| 00 | 03140 | 0.0880 
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Table 34. ARMA(3,4) MODEL, METHOD 3, POLE-ZERO COMPAR- 
ISON, REALIZATION 1. 


Zeros 
-0.4533 - 0.5691; 










b. Notse Realization 2 
Using the second noise realization, the coefficients converge after four it- 


erations. Fables 35 and 36 present the results. 


Table 35. ARMA(3.4) MODEL, METHOD 3, 
COEFFICIENTS COMPARISON, 
REALIZATION 2. 
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Table 36. ARMA(3,4) MODEL, METHOD 3, POLE-ZERO COMPAR- 
ISON, REALIZATION 2. 


Zeros 












The spectra using the true and the estimated network coefficients are plot- 
ted in Figure 17. 

When using method | the coefficients do not converge. Thev oscillate 
about two sets of values in an alternating fashion, hence the average of the two sets 1s 
used as the estimate of the coefficients in the comparison. The parameters of the first 
realization converge after nine iterations using method 2, and after six iterations using 
method 3. The parameters of the second realization converge after four iterations using 
method 3. Method 3 provides the best results in terms of the coefficients and pole-zero 
locations for both noise realizations. The coefficients converge also earlier when using 
method 3 compared to the other two methods. Method 2 performs not as well but gives 
also reasonable results. 

The spectrum due to the poles of the ARMA(3,4) models using the esti- 
mated coefficients closely resembles the original spectrum. The AR part coefficient es- 


timates are more accurate than the MA part coefficient estimates. 


D. THE ARMA(3,4) MODEL WITH OBSERVATION NOISE 

Observation noise is added to the ARMA(3,4) model and resulting noisy sequence 
is processed via the three different methods. The observation noise is independent of the 
driving noise. It is white Gaussian noise with zero mean and unit variance. The signal 
to noise ratio 1s about 15 dB. 


This model is illustrated in Figure 18. 
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Figure 16. The Spectra of ARNIA(3,4) Models, Method 2 
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Figure 17. = The Spectra of ARMA(3,4) Models, Method 3 
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Observation Noise 


Driving Noise 


E(n) 


ARMA(3,4) 





Figure 18. The ARMA(3,4) Model with Observation Noise 


1 METHOD 1 
The coefficients do not converge in nine iterations. Because of oscillations 
about two values, the average of the two values is used as an estimate of the coefficients 


in the comparison. Tables 37 and 38 present the results. 


Table 37. ARMA(34) MODEL WITH OB- 
SERVATION NOISE, METHOD 1, 
COEFFICIENTS COMPARISON. 


Estimated 
a, 2.0285 -0,4285 











i 9241 | 0.2559 
“1485 [0.0683 
dy 0.7225 0.8145 + 0.0920 









ae 
Toa Paras | #13813 
2.9370 
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Table 38. ARMA(3.4) MODEL WITH OBSERVATION NOISE, 
METHOD 1, POLE-ZERO COMPARISON 


Zeros 









The spectra using the true and the estimated network coefficients are plotted 
in Figure si) 
2, METHOD 2 
The coefficients do not converge in nine iterations. Because of oscillations 
about two values, the average of the two values is used as an estimate of the coefficients 


in the comparison. Tables 39 and 40 present the results. 


Table 39. ARNIA(3,4) MODEL WITH OB- 
SERVATION NOISE, METHOD 2, 
COEFFICIENTS COMPARISON. 


| 0.7235, 1.4933 | +0708 
61.00 | 8.5983 + 7.5983 
Tf 080 [3.434 3.0234 
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Table 40. ARMA(3,4) MODEL WITH OBSERVATION NOISE, 
METHOD 2, POLE-ZERO COMPARISON 


Poles- True Estimated Distance Radial Diff. 
Zeros 









The spectra using the true and the estimated network coefficients are plotted 
merrcure 20. 
3. METHOD 3 
The coefficients converge after four iterations. Tables 41 and 42 present the 


results. 


Table 41. ARMA(3.4) MODEL WITH OB- 
SERVATION NOISE, METHOD 3, 
COEFFICIENTS COMPARISON. 


a | 1.60 | 1.994] 40.1006 
a | 0.7235) 0.7353] 100128 
Tb 1.00 [3898 | 0.3898 
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Table 42. ARMA(3,4) MODEL WITH OBSERVATION NOISE, 
METHOD 3, POLE-ZERO COMPARISON 


Poles- True Estimated Distance Radial Diff. 
Zeros 









The spectra using the true and the estimated network coefficients are plotted 
in Figure 21. 

In summary, the estimated coefficients do not converge using method | or 2 
but they do converge after four iterations using method 3. The spectra using method | 
and 2 are poor. But the method 3 gives results close to the actual ones. The spectrum 
estimated with this method follows the original pattern except for the zero location. 


This 1s believed to be partially due to the imprecise MA coefficient estimation procedure. 
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ARMA(3,4) Model with Obser. Noise, Method 1 


——— Using true coeff. nlp 


see iUsing aver. of the est. coeff. 
s. ‘ ; iraundavarcoesaroua Tiers cares ie ' aisapea 6 


-o_~, 
Q 
= 
v 
a] 
es 
= 
=) 
bp 
0 
= 


Sainple Number 





Figure 19. The Spectra of ARNIA(3.4) Models with Obser. Noise, Method 1 
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ARMA(3,4) Model with Obser. Noise, Method 2 
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Figure 20. The Spectra of ARMA(3,4) Models with Obser. Noise, Method 2 
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ARMA(3,4) Model with Obser. Noise, Method 3 
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Figure 21. The Spectra of ARNIA(3,4) Models with Obser. Noise, Method 3 
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V. CONCLUSIONS AND RECOMMENDATIONS 


In this thesis, an iterative technique to estimate the coefficients of ARMA models 
driven by a random input is presented. The AR parameters are estimated initially, then 
MA parameters are estimated assuming the AR parameters have been compensated for. 
In an iterative fashion MA and AR contributions are removed from the original data 
allowing improved AR and MA coefficient estimates. Three different AR estimation 
methods are experimentally explored. The third method provides the best result in terms 
of the true and estimated coefficients and pole-zero locations. This method uses the 
pseudoinverse with correlation function values starting at the zero lag after removing the 
MA influence. For models of real time data which have an odd number of poles, one 
should address the issue of the DC component (i.e., remove it). Simulation results for 
an odd ordered pole model is presented in Appendix C. 

Also, by examining the spectra of the models, we can sav estimation of poles is 
obtained more accurately than the estimation of zeros. The Cholesky factorization is 
used for the estimation of the MA part coefficients but it is an approximate solution. 
For that reason, we do not expect superior results for the MA estimation part. 

Further research should concentrate on improving the MA part coefficients esti- 


mation. 
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APPENDIX A. COMPUTER PROGRAM TO ESTIMATE 
COEFFICIENTS 


WR ee ee ee ee 


% This program written by * 
% Gurhan Kayahan for IBM PC/AT * 


% Matlab package. % 
URE RE ee Te NTE Te Ve Fe Te Te ENE EVER TEN NE TENE TT NC TE 
PETE TTT TE TET TET TE TEE TEE TE TE ETE ELE TE TERETE TE TOTES TER TE 
% GENERATION o ARMA TIME SERIES * 
yas owen ute ton! onl ss esles ate seve ves Ye ste steste slew mt slew aio nfe steals sles aie stesles's sloule sleute ste seve ate 
t=sqrt(-1); 


rootsl=[ -—0.25+0.57 -—0.25-—-0.512];change zeros for each model 
roots2= [0.6+0.47 0.6-—0.47 0.8];change poles for each model 
b=poly(roots1); find true coefficients of MA part 
a=poly(roots2); find true coefficients of AR part 

a=real(a); 

t=1: 1000; 

y=t*0; 

rand( ‘normal'); generate random signal 

stddev=sqrt(1); 

rnum=stdev*rand(t); 

yi=filter(b,a,rnum); filter random signal through ARMA model 

x= [Lones(1,999)]; generate step signal 

y2=filter(b,a,x);find step response of ARMA model to discard transient response 
y=y1(52: 252)5 


Yer s\s s's sjeJe s'es'e Yesles'ss'e ses s'e ste slesie ste ateslente stes's steste sles'e slexts s'e sles'evlesie se sieve se ste siesie verx'e vee s'e steve se s'e toute sioviacte 
e ae AR COEFFICIENTS USING eee ees APPROACH % 
Brewers Vevesesevedotesevesesedesedetesededesetedtesedeesetededssesededeseseseeleteisieleteiedesleeieiedelesede 


Mes corr(y): 
iri 201: 399); 
r=flipy(rl); 
row=r( 194: 196); 
col=flipy(r( 192: 194)); 
a=toeplitz(col, row); 
Beelapy( -r( 191: 193)); 


c=b'; 

x=b*inv(a); square matrix inverse 

al=x'; 

Brereven sleste slew ao sie to slooten' steve sje sfc s'e Sexes! eveseves eves’. fo whe au! eute ste stes ate ste ve ves sleafevleslesies) ve sje sive ate lesles'c sfeateales slacieatacleute sesiesk slestesle 
% CALCULATE MA COEFFICIENTS USING CHOLESKY FACTORIZATION * 
Yu le ule sJes'e verte slesies slesies'es' eve ve os esfes evlon Se verte slevesie ve w‘e s'. eve sles ate ve ale viene ston ed ve ste ate ate le ele levee ents leslesteslestes esleve sles! lenteule site lente ate sto 
c= Ue 


n= [1 4l(1) @a1(2) 41(3)); 

@—11lter(n,d,y); 

Z1=z(4: 201); 

mex CcOrr( z1); 

mz=tZ \( 1985395); 

i=1; 

while i<199 
rz3(i)=rz(i)/198; 
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i=itl; 

end 

EZZ=L Zo le sake 

rzt=toeplitz(rz2); 

cho=chol(rzt); 
bl=c€ho( lao 
Wrere ved eds sete eee desededesesedededcadedevdeseseseacaedesdede Teves aes vededevese ve veve 
% ROUTINE FOR ESTIMATION OF AR COEFFICIENTS * 
YTV TEV NTE Te Ve Tee ve ge He ye Be Ve Te Ve Fe Te TENE Te Fe WEE VEE HE Vee Tee TENE Me Te Ve Vee TENE HE 
1=2; 
while i<ll 

d= [b4{7~—1,1) bl = 1,2) Sr ay 

ea 

Zar=tilter(n, day. 

rar=xcorr(zar); 

rl=rar(201: 399); 

r=f lipyGa: 


row=r( 194: 196); change correlation lags for each method 


col=flipy(r( 192: 194)); change correlation lags for each method 


g=toeplitz(col,row); 


c= flipy(- r(191:193)); change correlation lags for each model 


ia : 

x=f*inv(g); square matrix inverse (replace with 
VCs eo: 
al=.[al). aca an: 


ie wede pti se bore forts oe ich ein oe ee aise ae eS ee Fo SESS SESS LSS ais 
‘he pen ee eben OF MA ee ee ¥ 
Pe TeT ee TVET TER TE TE ETE TET ETE TE ETE TER TE TE TCDS ETE NE TELE TERETE ETE ERE ES CTEM 
o C1: 


n= £1 @(735)) 62, 2iaiz 5 3a 

Z—fLileer(n.d, ya: 

Z1=z(4: 201); 

EZ N= MeCOrr Zt) 

EoaEZlt 198; 29>) 

k=1 

while k<199 
rz3(k)=rz(k)/198; 
k=k+1; 

end 

EZO=TZ5 Glo 

rzt=toeplitz(rz2); 

cho=chol(rzt); 

Ae ec 
bl= [bY 6b Cine) 4. 


e(i-1)=sum( (al(i,: yeal(iel,:)). a 2 )/3+sum( (b1(i,: 


if e(i-1) < 0.0001 
al=al; 
bl=bl; 

else 
i otaelle 

end 

end 


x=pinv(g) for pseudo inverse) 


)-b1(i-1,: )) eee ee 


APPENDIX B. COMPUTER PROGRAM TO CALCULATE TRUE AND 
ESTIMATED SPECTRA 


WAM TEEN Te He Te TEN He Fe Ve Fe Te Te TE eK Fe FE Te Me Fe HE TENE TENS He HEE ETE 
° CALCULATION OF TRUE SPECTRUM 


Www Ce ek oie bie ok ae ek ir te ok oe elie ie tee ee ie ee ok ee ok et 
a= [1 -—1.20 0.52]; change true poles for each models 
b= [1 0.50 0.3125]; change true zeros for each models 
al= [é@ zeros({1,123)]; zero padding 
bi= (5 zeros(1,124)]; zero padding 
fal=fft(al); fast fourier transform of AR part 
fbl=fft(bl); fast fourier transform of MA part 
mfal=abs(fal). a 2; 
mfbl=abs(fbl). a 2; 
c=mfbl. /mfal; 
a=max(c); 
i=1; 
while i<129 
c(ij=c(i)/a; 
i=itl; 
end 
spec=10*logl0(c); 
VETTE TNT TE BRERA LE TC TE TERE TENET ETE EDL NETTIE TE 


%o ee ON OF iota Learn 


aes wevevevevevesedescsescdevesleseslesleveses eslesteveseseveslededlcslestesiccese 


a= 4 —1.1531 0.5013]; Eiicee Lee poles for each models 

b= (1.0174 0.4732 0.1202]; change estimated zeros for each models 

al= [a zeros(1,123)]; 

bl= (b zeros(1,124)}; 

fal=fft(al); 

Pel=tit(bl); 

mfal=abs(fal). a 

mfbl=abs(fbl). a 

c=mfbl. /mfal; 

a=max(c); 

i=1; 

while i<129 
e(ij=c(i)/a; 
i=it+1; 

end 

spec7= meee Cc); 


ws ale se s's Je sie ale ate se ate s'e Neste ves s'e ate owe wha aton' awe eveslesles! sledlestests ve se ste ve eseslesleslesics eskele Je we ese stevie 


% PLOTTING OF TRUE AND ESTIMATED SPECTRUM 

Wet stevlenteclentecleclovlolontes teslecles enlowlentacioclectoatsctouio stervlevievtes's steslevesleslesles' levevecles slonfontonte 

t=0: 127; 

t=t'; 

plot(t(1: 65) ,spec7(1:65),'*' ,t(1: 65),spec(1: 65)) 
fede (' ARMAC2 , Zz) Model, echoed 3, Realization 1') 
xlabel('Sample Number’ ) 

ylabel( 'Magnitude(dB)' ) 

grid 

text(25,-5,'---- Using true coefficients’) 


b) 


2 
a 
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text(25,-9, *%%" After 4 iterations’ ) 
delete specl.met 
meta specl 


APPENDIX CC. SIMULATION RESULTS FOR AN ARMA(Q,3) MODEL 


For models with an odd number of poles, one of the poles must be located on the 


real axis in the z-domain to generate real time data. 


The pole-zero locations for this model are illustrated in Figure 22. 





Figure 22. The ARNMA(2,3) Model Pole-zero Locations 


1. METHOD 1 
a. Notse Realization | 
The coefficients do not converge in nine iterations. Because of oscillations 
about two values, the average of the two values 1s used as an estimate of the coefficients 


in the comparison. Tables 43 and 44 present the results. 
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Table 43. ARMA(2,3) MODEL, METHOD 1, 
COEFFICIENTS COMPARISON, 
REALIZATION 1. 


Ta | 200 | 1.9995 | +0.0575 
| 1.00 | 1.0562 | +0.0562 
| 0.50 | osisa | 00184 














Table 44. ARMA(2,3) MODEL, METHOD 1, POLE-ZERO COMPAR- 
ISON. REALIZATION I. 


Poles- True Estimated Distance Radial Diff. 
Zeros 
0.6 + 0.4; | 0.6526 + 0.3809; 0.0596 
ae 0.6526 - 0.3809} 0.0596 


0.87 0.0000 
|, | -0.25 + 0.5) | -0.2454 + 0.3145; 0.1989 
02454-03145 0.1989 





b. Noise Realization 2 
Using a different noise realization, the coefficients converge after seven it- 


erations. Tables 45 and 46 present the results. 
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Table 45. ARMA(2,3) MODEL, METHOD 1, 
COEFFICIENTS COMPARISON, 
REALIZATION 2. 


a | 200 | 2.8571 | 0857 
3.7525_| 1.3725 


[1.00 [sas | +0348 
T4050 [oun | 03709 





Table 46. ARMA(2,.3) MODEL, METHOD 1, POLE-ZERO COMPAR- 
ISON, REALIZATION 2. 


Poles- True Estimated Distance Radial Diff. 
Zeros 


0.6 + 0.4) 0.8301 + 0.2765] 0.2611 0.2664 
0.6 - 0.4) 0.8301 - 0.2765} 0.2611 0.2664 


Tn fos .ises 0.3968 | 0.0000 
1071 
Te, [025-055 [-0.5266 





2, METHOD 2 
a. Noise Realization I 
The coefficients do not converge in nine iterations. Because of oscillations 
about two values, the average of the two values is used as an estimate of the coefficients 


in the comparison. Tables 47 and 48 present the results. 
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Table 47. ARMA(2,3) MODEL, METHOD 2, 
COEFFICIENTS COMPARISON, 
REALIZATION 1. 


0.3466 
97a | 0.507% 


[4.00 | 12677 | +0267 
4, | 050 | oszma | 403074 





Table 48. ARMA(2,3) MODEL, METHOD 2, POLE-ZERO COMPAR- 
ISON, REALIZATION 1. 


Zeros 

pm | 06-04) | 0.5267 - 0.2515) 
pos Por 0999 | 
Ts, [0.25 + 05] _|-0.3263 + OAR 
Te Ps 05j 0263-04; [016 | 0807 
















b. Notse Realization 2 
Using a second noise realization, the coefficients do not converge. Because 
of oscillations about two values, the average of the two values 1s used as an estimate of 


the coefficients in the comparison. Tables 49 and 50 present the results. 
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Table 49. ARMA(2,3) MODEL, METHOD 2, 
COEFFICIENTS COMPARISON, 
REALIZATION 2. 


“11003 | +0.8995 
‘oanaa | -1593: 


Ts, | 050 | 13061 | 0.8084 





Table 50. ARMA(2.3) MODEL, METHOD 2, POLE-ZERO COMPAR- 
ISON, REALIZATION 2. 


Poles- True Estimated Distance Radial Diff. 
Zeros 


0.6 + 0.4) | 0.8145 + 0.2883} 0.2418 0.2478 
0.8145 - 0.2883} 0.2418 


ep 
0.25 + 0.5 0.1919 | 0.3444 
Cayo) Oot? 0.3444 





3. METHOD 3 
a.  Notse Realization I 
The coefficients converge after five iterations. Tables 51 and 52 present the 


results. 
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Table 51. ARMA(2,3) MODEL, METHOD 3, 
COEFFICIENTS COMPARISON, 
REALIZATION 1. 


a | 200 | 1.4057 | +0808 
[a | oars | 0.0890 | +0.5050_ 


|b | 001.8594 | +0.8594 | 
| | 0.50 | 14885 | -0.9885_ 





Table 32. ARMA(2,.3) MODEL, METHOD 3, POLE-ZERO COMPAR- 
ISON, REALIZATION 1. 


Zeros 










b. Noise Realization 2 
Using a second noise realization, the coefficients do not converge. Because 
of oscillations about two values, the average of the two values is used as an estimate of 


the coefficients in the comparison. Tables 53 and 54 present the results. 
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Table 53. ARMA(2.3) MODEL, METHOD 3, 
COEFFICIENTS COMPARISON, 
REALIZATION 2. 


a [200 [1.7973 | +0200 
9.9959 | 0.8841 


Tf 1.00 [1.2530 | #02540 
Tf 03s] ons [0.2974 





Table 54. ARMA(2,3) MODEL, METHOD 3, POLE-ZERO COMPAR- 
ISON, REALIZATION 2. 


Poles- True Estimated Distance Radial Diff. 
Zeros 
0.6 + 0.43 0.1263 + 0.7679} 0.5997 0.2564 
0.6 - 0.4; 0.1263 - 0.7679} 0.5997 


1546 0.7446 | 0.0000 
-0.4645 0.5440 1.1071 
-0.0259 0.5479 1.1071 





Simulation results have shown that for the ARMA(2,3) model with one 
pole on positive real axis, realization 2 gives poor results. This is because the DC com- 
ponent for this realization is - 0.9017 while it is -0.0580 for realization |. For parameter 
estimations the DC component should be removed. Therefore the data is filtered effec- 
tively removing the pole on the positive real axis (i.e., DC component is removed and 
the ARMA(2,3) model becomes ARMA(2,2) model). The results of this ARMA(2.,2) 
model was presented in chapter IV-B. The ARMA(2,3) and ARMA(2,2) true and esti- 
mated spectra using method 1 are presented in Figure 23. Method 2 results are pre- 
sented in Figure 24 and method 3 results are presented in Figure 25. The simulation 


shows that the method 3 gives good results as long as the DC component is small. 


| 


ARMA(2.9) Model, Method f, Neallzntlosn § 


reeeensyy aon 
e 


— t ‘{felent ) 
~ oS eh alas okie 20S : wwest fisiug Irie cocfficients 
e ' 


sees Using avernge of the eat coefl. 


AMMA(Z,2) Model ft. Method 1, Neaffzation § 


0 *eg-snee -—— + 
e 


+ 
~ 


cece Afier 3 Hernatlons 
' 


Magnitude( dB) 


SO ee ee 


no 00 


Sample Numer Sample Number 


ARMA(2.3} Modc!, Method t, Realizallan 2 


+-----———- 


ARMA(2,2) Mode! 1, Method f, Realization 2 


wee fattg fror corffleients 


et gee. 


After)? ieralicns mame Uaing trie cocificlénts 


cece Afice 2 tterallons 


Magnitude(dB) 


ue 
40 


Sampte Number Sampte Number 





Figure 23. The Spectra of ARMA(2,3) and ARM A(2,2) Model-B, Method 1 
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Figure 24. The Spectra of ARMA(2,3) and ARMA(2,2) Model-B, Method 2 
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